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IV 


1.  SUMMARY 


Discriminants  based  on  the  differences  in  the  seismic  source  spectrum  between 
earthquakes  and  explosions  often  work  best  in  the  highest  recordable  frequency  band 
(e.g.,  Walter  et  al.,  1995;  Allmann  et  al.,  2008).  In  this  high  frequency  band,  (>  2-4  Hz) 
scattering  of  elastic  energy  by  small-scale  heterogeneity  (less  than  a  wavelength)  can 
equilibrate  energy  on  all  components  of  motion  and  stabilize  the  behavior  of  yield 
measures  based  on  the  Lg  wave  trapped  in  the  Earth’s  crust.  Lateral  variation  of  both 
small-scale  heterogeneity  and  larger-scale,  resolvable  structure  (greater  than  a 
wavelength)  can  control  the  blockage  of  Lg  and  the  amplitude  of  other  regional  seismic 
waves  (Cao  and  Muirhead,  1993;  Sens-Schonfelder  et  al.,  2009).  This  project  models  the 
combined  effects  of  the  large-scale  (deterministic)  and  the  small  scale  (statistical) 
structure  to  invert  for  improved  structural  models  and  to  evaluate  the  performance  of 
yield  estimators  and  discriminants  at  selected  seismic  stations  in  Eurasia.  These  tasks  are 
approached  by  synthesizing  seismograms  using  a  radiative  transport  technique  to  predict 
the  high  frequency  coda  (>5  Hz)  of  regional  seismic  phases  at  stations  having  known 
large-scale  three-dimensional  structure,  combined  with  experiments  to  estimate  the 
effects  of  multiple-scattering  from  unknown  small-scale  structure.  We  describe  a 
radiative  transport  code  and  report  the  results  of  its  use  in  modeling  regional  wave 
propagation,  including  tests  of  sensitivity  of  the  shapes  of  high  frequency  seismic  coda 
shapes  to  structural  and  source  parameters  in  the  Lop  Nor,  China  region. 

2.  INTRODUCTION 

This  section  motivates  the  needs  of  CTBT  monitoring  for  modeling  high-frequency 
seismic  waveforms  in  3-D  structures  at  ranges  much  longer  than  100  wavelengths. 

Section  3  describes  how  a  radiative  transport  algorithm  can  fulfill  this  need  and  outlines 
our  development  and  implementation  of  such  an  algorithm  with  an  open-source  computer 
code.  Section  4  presents  results  of  experiments  that  vary  deterministic  and  statistical 
parameters  of  Earth  models  as  well  as  source  radiation  patterns,  depth,  and  tectonic 
release.  Concluding  remarks  are  presented  in  Section  5. 

2.1  The  importance  of  very  high  frequency  coda  modeling 

Differences  in  source  properties  between  a  theoretically  ideal  earthquake  and 
explosion  are  well  understood:  distributed  slip  on  a  fault  plane  versus  an  approximate 
step-function  in  pressure  applied  to  a  spherical  cavity  deep  enough  for  explosion 
containment  (Brune,  1970;  Mueller  and  Murphy,  1971).  These  simple  theories  do  not 
incorporate  dynamically  triggered  fault  motion,  complex  surfaces  of  equivalent  explosion 
volumes,  and  explosive  tectonic  release.  Nevertheless  they  are  good  enough  to  accurately 
predict  that  an  earthquake  with  the  same  scalar  moment  of  will  have  much  higher  comer 
frequency  in  its  displacement  spectrum.  The  success  of  common  seismic  discriminants 
based  on  amplitude  ratios  in  different  group  velocity  windows  (e.g.,  Pg/Lg)  depend  as 
much  on  this  spectral  difference  as  they  do  on  differences  in  source  depth  and  P  and  S 
radiation  patterns  (Walter  et  al,  2009). 
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High  frequency  Lg  coda  has  also  been  shown  to  be  useful  for  yield  estimation  (e.g., 
Mayeda,  1993).  As  frequency  increases,  however,  it  becomes  increasingly  difficult  to 
obtain  perfect  knowledge  of  small-scale  structure  needed  to  predict  a  seismogram  of  an 
earthquake  or  an  explosion  wiggle  for  wiggle.  Yet  even  without  this  knowledge,  simple 
statistical  descriptions  of  medium  properties  can  be  used  to  predict  the  shape  of  high 
frequency  Lg  coda. 

2.2  Limitations  of  numerical  modeling:  high  frequency  and  long  range 

Given  the  importance  of  high  frequencies  to  seismic  monitoring,  it  would  be  desirable 
to  model  high  frequency  regional  seismograms  with  enough  computational  efficiency  that 
many  structural  and  source  models  could  be  routinely  tested  and  compared  against 
observations  within  several  hours  or  less.  Unfortunately,  even  with  the  advance  of 
parallel  computation  and  increased  access  to  large  clusters  of  processors,  numerical 
modeling  in  three-dimensional  earth  models  rarely  exceeds  ranges  of  several  hundred 
wavelengths  for  computing  times  less  than  1  day.  At  5  Hz  this  range  is  typically  less  than 
100-200  km,  which  is  valuable  to  predictions  of  earthquake  strong  ground  motion,  but 
more  limited  in  value  to  nuclear  monitoring  where  IMS  stations  may  be  1000  km  or 
further  from  regions  of  interest  in  Eurasia. 

2.3  Importance  of  large-scale  structure  and  limitations  of  path 
calibrations 

Given  the  apparent  success  of  discriminants  and  yield  estimators  from  high  frequency 
coda  with  minimal  knowledge  of  earth  structure  and  source  descriptions,  one  may  believe 
that  monitoring  is  purely  a  problem  of  calibrating  coda  measurements  and  amplitude 
ratios  along  paths  to  each  monitoring  station,  applying  spatial  interpolators  (e.g,  kriging, 
Fan  et  al.,  2002)  and  a  model  of  propagation  in  a  plane-layered  medium  (e.g.,  MDAC, 
Taylor  et  al.,  2002).  Many  examples,  however,  have  been  collected  in  which  sharp 
gradients  in  the  Moho  depth,  mountain  roots,  and  deep  sedimentary  basins  have  strongly 
affected  azimuthal  variations  in  the  character  and  detection  of  regional  seismic  phases 
(e.g.,  Ruzaikin  et  al.,  1971),  making  simple  propagation  models  and  path  interpolation 
unreliable  tools  for  predicting  waveform  behavior. 

The  effects  of  larger-scale  (greater  than  a  wavelength)  structure  will  hence  limit  the 
portability  of  discriminants  and  yield  estimators  and  limit  the  effectiveness  of  their  path 
interpolations  and  calibrations.  In  addition  to  these  operational  concerns,  a  number  of 
fundamental  problems  involving  the  interaction  of  large  with  smaller-scale  (wavelength 
and  smaller)  structure  have  also  not  been  well  investigated: 

•  Can  large-scale  structures  always  block  the  highest  detectable  frequency  of 
some  regional  phases? 

•  Is  there  a  frequency  band  on  Earth  in  which  the  effects  on  regional 
seismograms  of  large-scale  structure  disappear? 

•  How  important  is  small-scale  structure  in  the  vicinity  of  large-scale  structure? 
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•  Is  anisotropy  of  the  scale-lengths  of  small-scale  structure  (horizontally  or 
vertically  stretched  heterogeneity)  important  to  regional  phase  propagation? 

By  combining  the  effects  of  large-scale  and  small-scale  structure  in  a  computationally 
efficient  3-D  method  of  seismogram  synthesis  based  on  the  radiative  transport  approach, 
our  project  seeks  to  determine  small-and  large-scale  3-D  structure  in  the  vicinity  of 
selected  monitoring  stations  by  matching  the  coda  of  very  high  frequency  (>  5  Hz) 
regional  phases.  The  results  of  this  modeling  will  be  used  to  optimize  coda-based  yield 
estimators  and  discriminants  as  a  function  of  azimuth  to  these  stations. 

3.  TECHNICAL  APPROACH 
3.1  Radiative  Transport 

Radiative  transport  seeks  to  model  the  diffusive  transport  of  energy  through  processes 
of  multiple  scattering  in  a  medium  having  random  fluctuations  in  physical  properties.  Its 
earliest  and  most  extensive  applications  are  to  the  propagation  of  electromagnetic  waves, 
especially  to  the  description  of  stellar  images  and  scintillation  (Chandrasekhar,  1960). 

Wu  (1985)  introduced  radiative  transport  to  seismology,  modeling  the  coda,  or  envelope 
of  energy,  starting  and  following  elastic  P  and  S  waves.  By  including  multiple-scattering, 
radiative  transport  improved  upon  the  single  scattering,  Rayleigh-Born,  model  of  seismic 
coda  introduced  by  Aki  and  Chouet  (1975).  An  early  significant  result  of  Wu’s  work  with 
radiative  transport  is  that  multiple  scattering  often  becomes  the  dominant  mechanism 
controlling  coda  shapes  at  frequencies  above  5  Hz  on  Earth,  explaining  a  spindle-shaped 
(lunar-like)  coda  observed  in  this  frequency  band.  Reviews  of  seismic  radiative  transport 
modeling  include  the  text  by  Sato  and  Fehler  (1998)  and  a  monograph  by  Margerin 
(2004).  Przybilla,  et  al.  (2009)  has  compared  the  method  against  finite-difference 
synthetics  in  2-D,  validating  body  wave  coda  predictions  of  the  radiative  transport 
method. 

Following  Shearer  and  Earle  (2008),  an  algorithm  to  implement  radiative 
transport  can  be  described  in  a  few  simple  steps: 

(1)  Choose  a  heterogeneity  power  spectrum  for  the  medium  (the  heterogeneity  power 
spectrum  is  the  Fourier  transform  of  the  spatial  autocorrelation  of  small-scale 
heterogeneity). 

(2)  From  the  heterogeneity  spectrum,  calculate  the  scattering  coefficient,  g,  which  is 
the  scattering  power  per  unit  volume. 

(3)  Average  g  for  P  and  S  waves  over  all  angles,  forming,  for  example,  g0  =  l/lp, 
where  path  lp  is  the  mean  free  path  of  a  P  wave. 
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(4)  Calculate  path  length  r,  assuming  r  is  an  exponentially  distributed  random  number 
having  a  mean  value  lp  for  P  waves  and  ls  for  S  waves  respectively,  i.e,  rp  =  -ln(x), 
where  x  is  a  random  number  between  0  and  1 . 

(5)  Choose  the  scattered  direction  from  sampling  a  probability  density  taken  from  the 
Bom  scattering  coefficients  of  the  heterogeneity  model  (Sato  and  Fehler,  1998) 
and  continue  computing  the  path  of  the  scattered  energy  packet.  Random  numbers 
sampling  a  scatter  radiation  pattern  are  used  to  determine  whether  the  scattering  is 
into  P  or  S  and  to  determine  S  polarization. 

Figure  1  shows  how  multiply  scattered  packets  of  elastic  energy  are  tracked  starting 
and  their  relation  to  a  statistical  realization  of  small-scale  heterogeneity. 


Itroei'jmr 

V 


{b)  won  Karma n  *  -  1  fl 


[c>  von  Karma n 
k  =  0.5  (Exponential) 


[dj  von  Knrirmn  k  «0/i 


Figure  1 .  Radiative  transport  algorithm  and  heterogeneity  spectra.  Top  left:  random  walk 
along  ray  paths  of packets  of  elastic  energy  scattered  by  heterogeneity  having  a 
wavenumber  spectrum  (top  right)  parameterized  by  a  scale  length  a  after  which  small 
spatial  scales  by  varying  powers  ofka  (ma)  where  k  =  2  jz/wa  velength.  Bottom:  example 
spatial  realizations  of  heterogeneity  spectra. 
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3.2  Advantages  of  Radiative  Transport 

The  principle  advantage  of  radiative  transport  is  the  ability  to  simulate  the  effects 
of  small-scale  structure  without  the  need  for  model  meshes  that  are  dense  enough  to 
explicitly  describe  that  structure.  Instead,  the  model  mesh  describes  only  the  large-scale, 
long-wavelength  equivalent,  background  medium,  against  which  small-scale 
heterogeneities  are  assumed  as  a  random  perturbation  field.  The  small-scale 
heterogeneities  are  described  in  the  aggregate  via  a  small  number  of  statistical 
parameters,  rather  than  describing  the  perturbation  field  explicitly  in  fine  detail.  This 
represents  a  dramatic  reduction  in  the  storage  requirements  needed  to  describe  the  Earth 
model.  Additionally,  since  radiative  transport  is  not  a  numerical  wave-equation  solver,  it 
also  avoids  the  requirement  for  a  simulation  using  a  grid  dense  enough  to  capture  the 
temporal  and  spatial  derivatives  in  the  full  equations  of  motion.  Thus,  models  are 
described  with  what  we  call  a  “model  mesh”  rather  than  a  “simulation  grid,”  where  the 
model  features  determine  the  necessary  mesh  density  we  wish  to  capture,  rather  than  the 
wavelengths  we  wish  to  simulate.  In  contrast,  high-order  numerical  finite  difference 
methods  might  require  as  many  as  10  to  50  grid  points  per  wavelength.  This  would 
quickly  become  a  limiting  factor  on  the  spatial  extent  of  models,  making  unfeasible 
routine  modeling  of  seismic  waves  that  exceed  several  100  wavelengths  of  3-D 
propagation 

This  sparse-mesh  advantage  depends  on  the  ability  to  categorically  separate 
structure  into  large  and  small  scales.  In  general,  the  distinction  is  made  by  comparison  of 
the  feature  size  to  the  wavelength  being  simulated.  Large  structure  is  that  for  which  ka 
»  1,  where  a  is  a  typical  feature  size  and  k  is  the  wavenumber  corresponding  to  the 
chosen  wavelength,  and  small  scale  is,  for  radiative  transport  purposes,  structure  for 
which  ka  ~  1,  which  is  the  regime  in  which  scattering  is  most  pronounced.  Radiative 
transport  simulations  are  usually  carried  out  at  a  single  frequency,  and  this  choice  of 
frequency  establishes  the  reference  size  for  this  distinction.  Our  frequency  range  of 
interest  is  the  1  Hz  to  10  Hz  band,  which  at  typical  seismic  velocities  of  3  to  8  km/s, 
establishes  “small  scale”  in  the  approximate  range  of  0.3  km  to  8.0  km  or  smaller,  and 
means  that  node  spacing  on  our  model  meshes  can  be  much  larger  than  that. 

For  practical  modeling  purposes,  large-scale  structure  is  that  which  comes  from 
published  velocity  models  of  the  Earth,  in  which  seismic  wave  velocities  have  been 
imaged  from  receiver  functions,  travel-times,  and  waveforms  inversions.  Typically,  these 
images  are  determined  from  longer  wavelength  signals  than  the  high-frequency  signals 
that  are  of  interest  regional  nuclear  monitoring.  These  inversions  image  the  long- 
wavelength  equivalent  or  composite  medium  average  seismic  velocities.  These  large- 
scale  velocity  maps  then  form  the  background  structure  against  which  small-scale 
heterogeneities  are  assumed  as  a  perturbation  field,  for  use  in  generating  higher- 
frequency  synthetics.  The  advantage  of  radiative  transport  is  the  ability  to  describe  this 
perturbation  field  statistically  and  to  simulate  it  stochastically,  rather  than 
deterministically,  even  as  the  large-scale  structure  is  handled  by  deterministic  ray  tracing 
in  the  background,  long-wavelength  equivalent,  model. 
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The  practical  mesh  fineness  required  to  describe  deterministic  structure  is 
application-specific,  and  can  vary  quite  a  bit.  Node  spacings  of  10’s  to  100’s  of 
kilometers  for  lateral  variations,  and  five  to  20  kilometers  for  vertical  variation,  could  be 
considered  well  within  the  range  of  typical.  Scales  as  small  as  1  km,  however,  might  be 
needed  to  describe  lateral  variations  in  important  sharp,  well-resolved,  crustal 
discontinuities,  such  as  the  Moho.  A  typical  scale  length  for  the  material  heterogeneities 
that  become  the  stochastic  perturbation  field  might  easily  span  0.01  to  10.0  kilometers.  If 
features  of  these,  unresolved,  smaller  scale  lengths  were  to  be  described  explicitly,  the 
model  mesh  storage  requirements  would  increase  by  many  orders  of  magnitude. 

3.3  Construction  of  Models 

Earth  models  consist  of  a  background  deterministic  component,  consisting  of 
first-order  discontinuities  with  topography  separating  layers  having  a  3-D  variation  of 
and  velocities  and  densities  interpolated  over  a  course  grid,  and  a  perturbed  statistical 
component,  specified  by  a  heterogeneity  spectrum  having  parameters  specifying  the 
power  of  velocity  and  density  fluctuation  and  the  shape  of  this  spectrum  as  a  function  of 
wavenumber  (Figures  1  and  2). 

Frankel  and  Clayton  (1986)  demonstrated  how  earth  models  can  be  constructed 
that  have  a  specified  heterogeneity  spectrum.  The  procedure  consists  in  having 
perturbations  at  knot  points  driven  by  a  random  number  generator,  Fourier  transforming 
the  spatial  model  into  the  wavenumber  domain,  filtering  by  a  desired  wavenumber 
spectral  shape,  and  inverse  Fourier  transforming  back  to  space. 

Since  Frankel  and  Clayton’s  original  work,  advances  have  occurred  in  statistical 
characterization  and  understanding  of  small-scale  geologic  heterogeneity.  Goff  et  al. 
(1994)  described  a  procedure  by  which  models  having  the  spatial  statistics  of 
polycrystalline  or  multi-modal  assemblages  of  rocks  can  be  generated.  Work  beginning 
with  Fevander  et  al.  (1994)  and  Pullammanappallil  et  al.  (1997)  makes  it  possible  to 
formulate  statistical  models  for  common  sedimentary  and  metamorphic  formations. 

Model  statistics  can  be  modified  as  needed  to  reproduce  the  three-component  behavior  of 
the  coda  of  regional  seismic  phases.  Strong  emphasis  is  placed  on  correctly 
characterizing  the  statistics  of  small-scale  heterogeneities  in  the  upper,  highly 
heterogeneous,  5  km  of  the  earth.  Small-scale  statistics  affect  the  partitioning  of  P  and  S 
energy  on  the  three-components  of  motion  by  the  effects  of  scattering  near  the  source  and 
receiver  and  hence  the  performance  of  discriminants  tuned  to  differences  in  sources 
occurring  at  the  depths  of  contained  nuclear  tests. 

Lateral  variations  in  crustal  thickness,  basin  depths,  mountain  roots,  and  lateral 
tectonic  transitions  significantly  affect  the  phases  used  for  discrimination  and  detection. 
One  example  is  the  study  by  Pedersen  et  al.  (1998),  who  explain  anomalous  Rayleigh  to 
Love  mode  conversion  from  Lop  Nor  explosions  by  changes  in  crustal  thickness  at  the 
boundary  of  the  Tarim  Basin  and  Tian  Shan  mountain  belt.  Moho  topography  and  basin 
thickness  can  also  strongly  affect  the  propagation  of  Lg  (Cormier  and  Anderson,  2004). 
The  scale  of  these  types  of  lateral  structural  variations  is  often  large  enough  to  be 
resolvable  by  local  and  regional  reflection  and  refraction  experiments,  gravity  and 
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magnetic  data,  regionalization  by  surficial  geology,  and  global  surface  wave  inversions, 
e.g.,  Figure  2.  Hence,  we  refer  to  these  types  of  structures  as  deterministic.  The  types  of 
data  used  to  infer  deterministic  structure  are  collected  at  widely  different  spatial  scales, 
presenting  a  challenge  to  the  parameterization  of  a  three-dimensional  model  appropriate 
for  a  region  surrounding  a  particular  seismic  station.  The  parameterization  should  be 
flexible  enough  to  be  specified  at  high  resolution  where  data  is  available  and  at  lower 
resolution  where  it  is  not.  Resolution  should  be  high  to  describe  features  important  to 
regional  wave  propagation,  such  as  Moho  and  basin  topography,  but  can  be  lower  near 
interfaces  having  smaller  velocity  contrasts  and  lower  with  increasing  depth  in  the 
mantle,  where  heterogeneity  power  decreases. 

DETERMINISTIC  STRUCTURE 

Examples: 

•  Changes  in  Moho  depth 


•  Lateral  variation  in  seismic  velocit 


STATISTICAL  STRUCTURE 

Example: 

•  fine-scale  deviations  of  seismic 
velocity,  due  to  material 
inhomogeneity,  small  cracks  and 
fissures,  etc.  Random  heterogeneity 
can  be  parameterized  by  scale- 
length  and  strength  parameters. 


Figure  2.  Examples  of  deterministic  and  statistical  structure.  Top:  examples  of  known 
(deterministic)  large-scale  heterogeneity  and  its  effects  on  S  and  Lg  waves  (multiple 
critically  reflected  S  waves  at  the  Moho).  Bottom:  example  of  a  realization  of  3-D  small- 
scale  (statistically  described)  heterogeneity. 
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3.4  Software  Package 

3.4.1  Implementation  and  relation  to  existing  codes 

We  have  developed  a  new  software  package  for  synthetic  seismogram  generation 
called  Radiative3D  —  A  code  for  radiative  transport  in  3D  Earth  models.  (Homepage: 
https://rainbow.phys.uconn.edu/geowiki/Radiative3D)  Radiative3D  is  a  comprehensive 
new  code  with  flexible  Earth  model  entry,  moment-tensor  based  source  modeling,  and 
receiver  and  array  modeling.  The  simulation  engine  includes  deterministic  ray  tracing  for 
longitudinal  and  shear  propagation  modes  in  linear  velocity  gradients,  stochastic 
scattering  controlled  by  heterogeneity  parameters,  reflection  and  transmission  handling  at 
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sharp  boundaries,  and  attenuation  modeling  from  internal  friction.  Visualization  and 
post-processing  tools  have  been  created  to  produce  seismic  envelopes,  travel-time  curves, 
and  three-dimensional  body-wave  visualizations  of  energy  propagation  from  Radiative3D 
output.  Radiative3D  will  be  made  available  to  the  scientific  community  as  a  free  and 
open-source  software  product.  Radiative3D  is  a  new  code  that  takes  major  inspiration 
from  two  prior  codes,  namely  Raytrace3d  (Menke,  2002)  and  PSPhonon  (Shearer  and 
Earle,  2004).  Raytrace3D  solves  ray  propagation  in  a  3-D  velocity  model,  and  PSPhonon 
performs  radiative  transport,  including  scattering,  in  a  1-D  layered-earth  model.  Our 
code  performs  both  tasks  in  a  full  3-D  model. 

In  our  implementation,  elastic  energy  is  considered  to  originate  from  a  source  event  at 
a  fixed  location  and  time.  A  quantum  of  that  energy  is  then  propagated  as  a  discrete 
packet  or  bundle,  which  we  refer  to  as  a  phonon  (in  analogy  to  the  particle  representation 
of  light  as  a  stream  of  photons).  A  phonon’s  path  through  an  Earth  model  is  determined 
by  a  combination  of  ray  theory  to  handle  the  deterministic  propagation  through  the 
composite  medium  background  structure  (large  scale  structure,  capturing  broad  variations 
in  elastic  properties),  and  by  scattering  theory,  which  is  the  stochastic  handling  of 
scattering  due  to  small-scale  heterogeneities,  assumed  as  perturbations  to  the  background 
structure.  A  scattering  event  is,  in  essence,  an  interruption  and  randomization  of  a 
particle’s  otherwise  deterministic  progress.  With  a  sufficient  number  of  phonons  emitted 
from  the  source  event,  a  picture  of  the  energy  transport  throughout  the  model  begins  to 
emerge.  Records  of  the  phonon  travel  paths  can  be  used  either  in  bulk  or  via  selection 
criteria  to  visualize  this  energy  transport.  In  bulk,  the  data  can  be  used  to  produce  movies 
or  still-frames  of  the  evolving  wavefront  throughout  the  three-dimensional  model.  Via 
selection  criteria,  just  those  phonons  that,  for  example,  interact  with  the  surface  within  a 
specified  gather  radius  of  a  hypothetical  seismometer  can  be  collected  and  used  to 
produce  seismic  waveforms  of  surface  movement  at  the  given  location.  With  a  sufficient 
quantity  of  virtual  seismometers,  travel-time  curves  can  also  be  produced. 


3.4.2  Capabilities  of  Radiative3D 

Source  modeling:  Radiative3D  uses  a  moment-tensor  representation  to  describe 
seismic  sources.  Moment  tensor  representation  is  a  truncation  at  quadrupole  order  of  the 
multipole  expansion  of  the  elastic  wave  field  at  large  distances  from  the  source,  and  thus 
provides  a  way  to  parameterize  complex  sources  as  point-source  equivalents.  Despite 
discarding  higher-order  multipoles,  moment  tensors  provide  good  representation  of  a 
wide  range  of  naturally  occurring  and  man-made  seismic  sources.  Both  idealized  shear 
dislocation  (e.g.  slip  on  a  fault  plane)  and  idealized  explosions  are  exactly  representable 
by  moment  tensors.  In  fact,  the  six  degrees  of  freedom  in  the  rank-2  symmetric  moment 
tensor  allow  for  the  specification  of  a  full  range  of  event  types  from  pure  isotropic 
(explosion/implosion),  to  double  couple  (shear-dislocation  earthquake),  and  various  other 
configurations  such  as  compensated  linear  vector  dipole  (CLVD),  etc.  Internally,  the 
moment  tensor  is  converted  into  probability  distributions  describing  the  likelihood  of 
phonon  emission  at  various  take-off  angles  from  the  event  source. 
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Scattering  shapes:  A  scattering  event  involves  an  interruption  of  the  phonon  ray 
path  at  a  scattering  location  and  a  re-radiation  of  the  phonon  at  a  deflected  angle.  The 
angular  dependence  of  scattering  is  determined  by  the  scattering  radiation  patterns 
described  in  Sato  and  Fehler  (1998,  chapter  4).  These  radiation  patterns  depend  on 
frequency  and  on  four  statistical  parameters  describing  the  heterogeneity,  as  well  as  on 
the  background  P  and  S  elastic  velocities.  The  four  statistical  parameters  include  the 
strength  e  of  the  perturbation  field,  the  scale  length  a  that  defines  a  low-pass  comer  in  the 
power  spectrum  of  length  scales,  a  parameter  k  that  determines  the  fall-off  rate  of  the 
power  spectrum,  and  a  parameter  v  that  establishes  a  ratio  between  velocity  and  density 
fluctuation  strengths.  Like  event  sources,  the  scattering  shapes  are  also  represented 
internally  as  probability  distributions  characterizing  the  relative  likelihood  of  various 
deflection  angles. 

Flexible  Earth  model  input:  Radiative3D  supports  Earth  models  based  on  a 
multitude  of  model  construction  metaphors,  allowing  both  very  simple  and  very  complex 
models  to  be  constructed  with  ease.  One  supported  model  type  is  the  layered-Earth 
model,  in  which  the  model  consists  of  consecutive  layers  at  increasing  depth.  Unlike 
traditional  layered  Earth  models,  however,  a  degree  of  lateral  variation  is  supportable  by 
the  ability  of  the  interfaces  between  layers  to  take  on  arbitrary  orientation.  Thus,  for 
example,  gradually  broadening  or  gradually  tapering  crust  structures  can  be  modeled.  In 
layered  models,  the  material  properties  (elastic  velocities,  density,  heterogeneity  spectra, 
and  intrinsic  attenuation)  are  uniform  throughout  the  layer.  Gradients  in  properties,  e.g. 
velocity  gradients,  must  be  simulated  in  stair-step  fashion.  An  Earth-coordinate 
subsystem  (ECS)  allows  model  parameters,  event  locations,  and  seismometer  locations  to 
be  specified  in  a  coordinate  system  of  the  user’s  choice  and  automatic  conversion  to 
layered  internal  model  space  will  occur.  The  ECS  also  allows  Earth-flattening 
transformations  to  be  applied  as  a  user  option,  to  simulate  Earth  curvature  in  layered 
models. 

For  more  sophisticated  topographical  or  laterally  varying  models,  Radiative3D  also 
supports  models  composed  of  tetrahedra,  allowing  full  three-dimensional  models  to  be 
specified.  In  tetrahedral  models,  properties  are  specified  on  the  comers  of  the  tetrahedra, 
and  are  linearly  interpolated  throughout  the  volume  of  the  tetrahedron.  Thus  velocity 
gradients  can  be  simulated  directly.  Additionally,  values  are  continuous  across 
tetrahedral  boundaries,  unless  the  boundary  represents  a  sharp  discontinuity,  which  is 
also  supported  by  the  software.  At  present,  the  tetrahedral  functionality  of  Radiative3D 
has  been  recently  completed  but  is  still  in  testing  phase.  Tools  to  help  the  construction  of 
complex  models  for  input  to  the  software  are  also  under  development. 

Signal  output:  Radiative3D  has  two  basic  output  modes.  The  first  is  a  play-by-play 
event  stream  that  reports  everything  that  happens  to  a  phonon  from  generation  at  the 
source  event  until  eventual  phonon-death  as  the  phonon  either  leaves  the  model 
boundaries,  or  is  abandoned  after  its  lifetime  exceeds  the  window  of  interest.  This  mode 
can  be  used  for  volumetric  visualization  of  energy  propagation  throughout  the  whole 
model.  The  second  output  stream  is  the  time-binned  signals  from  virtual  seismometers 
that  collect  and  report  on  the  energy  accumulated  at  fixed  points  on  the  model  surface. 
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These  can  be  used  to  generate  synthetic  seismograms  (as  energy  envelopes)  or  synthetic 
travel-time  curves. 

3.4.3  Visualizing  Output 


Radiative3D  produces  output  suitable  for  visualization  with  external  tools.  The 
output  capabilities  fall  into  two  major  categories:  (1)  event  reporting,  and  (2)  seismic 
energy  binning.  Event  reporting  means  reporting  the  progress  of  individual  phonons  in  a 
play-by-play  manner.  Examples  of  “events”  in  this  context  include:  generation, 
scattering,  crossing  a  model  boundary,  free-surface  or  discontinuity  interface  reflection 
and  transmission,  etc.  When  these  events  are  analyzed  in  post-processing,  detailed 
pictures  of  energy  propagation  can  be  constructed.  One  of  the  first  visualization  tools 
developed  was  a  GNU  Octave  script  to  produce  videos  illustrating  P  and  S  energy 
propagation  in  the  Earth  model. 

Figure  3  shows  energy  propagation  visualized  as  a  time  series  plot  of  phonon 
locations  in  a  model  cross  section.  Red  dots  represent  P  phonons,  blue  dots  represent  S 
phonons.  Reflection  and  refraction  occurs  at  interfaces  between  model  layers  at  which 
velocities  are  discontinuous.  Conversions  between  P  and  S  polarization  are  driven  by 
both  scattering  and  reflection/transmission  events. 


Earthquake  Time-Series: 


Explosion  Time-Series: 


Figure  3.  Time  snapshots  of  regional  seismic  wavefield  versus  range  and  depth. 
Energy  propagation  visualized  as  a  time  series  plot  of  phonon  locations  in  a  model 
cross  section.  Note  the  initial  absence  of  the  S  wavefronts  (blue)  in  the  explosion 
time-series  and  its  progressive  development  with  time  due  to  multiple  scattering  and 
interaction  with  the  free  surface. 
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3.4.4  Seismometer  Output 


The  other  form  of  output  which  can  be  analyzed  for  quantification  or  visualization 
is  binned  seismic  energy.  Radiative3D  allows  for  the  placement  of  virtual  seismometers 
along  any  surface  designated  as  a  collection  surface.  Whenever  a  phonon  interacts  with  a 
collection  surface  within  a  specified  gather  radius  of  a  seismometer,  the  energy  of  that 
phonon  is  collected  by  the  seismometer,  decomposed  into  three  cartesian  energy  axes, 
and  binned  into  time  windows.  At  the  end  of  simulation,  these  energy  bins  are  output  and 
can  be  interpreted  as  seismic  energy  traces,  suitable  for  making  envelope  plots.  Figure  4 
shows  two  seismic  envelope  plots  produced  by  the  Radiative3D  code.  The  traces 
represent  2.0  Hz  phonon  energy  arriving  at  a  virtual  seismometer  approximately  800  km 
from  a  hypothesized  earthquake  source  (left)  and  explosion  source  (right)  synthesized  in 
a  3-D  model  of  the  Lop  Nor  region  having  a  simple  deterministic  and  trial  statistical 
structure.  We  set  the  depth  of  the  explosion  event  to  2.0  km  below  the  surface,  and  the 
depth  of  the  earthquake  event  to  42.0  km  below  the  surface,  comparable  to  some  of  the 
larger  magnitude  earthquakes  to  recorded  in  the  Lop  Nor  region. 

Traces  from  multiple  seismometers  arranged  in  a  linear  array  can  be  combined  to 
make  travel-time  curves.  Figure  5  shows  two  travel-time  curves  produces  by 
Radiative3D.  The  curves  represent  2.0  Hz  phonon  energy  from  earthquake  and  explosion 
sources  in  a  Lop  Nor  model  with  complexity  near  the  Moho  interface  and  a  velocity  ramp 
just  underneath  it,  resulting  in  phase  arrivals  with  Pn  and  Sn  timing  characteristics. 
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Figure  4.  Example  coda  envelopes  for  earthquakes  and  explosions  synthesized  by 
radiative  transport.  Seismic  envelope  plots  produced  by  the  Radiative3D  code  for  a 
shallow-focus  earthquake  (left)  and  explosion  (right).  Traces  are  shown  for  components 
of  motion  NS  (Y),  EW  (X),  and  vertical  (Z)  and  wave  type  (P  or  S)  arriving  at  the 
receiver.  Shaded  area  under  P  and  S  curves  represents  phonon  capture  count  and  serves 
a  diagnostic  purpose. 
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Figure  5.  Example  coda  envelopes  displayed  as  travel  time  versus  distance  plots.  Travel 
time  curves  produced  by  the  RadiativeSD  code  for  a  shallow  focus  earthquake  (left)  and 
explosion  (right).  Image  density  indicates  total  energy  arriving  on  all  three  axes,  and 
makes  various  seismic  phases  easily  discernable.  Left:  Earthquake  source.  Right: 
Explosion  source. 


3.4.5  Features  Complete: 

Layered  multi-cell  Earth  models: 

While  full-3D  Earth  models  based  on  tetrahedral  model  cells  are  not  yet 
supported,  an  intermediate  cell  type  based  on  stacked  cylinders  is  supported  and 
functional  for  simulating  layered-Earth  models.  Partial-3D  functionality  is  available 
because  the  interfaces  between  adjacent  cylinders  can  be  given  an  arbitrary  inclination. 
This  means  that  first-order  lateral  variations,  such  as  a  sloped  Moho  can  be  modeled. 

Reflection/Transmission  scattering  at  discontinuity  interfaces: 

Free-surface  reflections,  as  well  as  sharp  velocity  discontinuities  (such  as  at  the 
Moho  layer),  require  special  treatment  of  ray  propagation.  Much  like  the  radiation 
patterns  of  scattering  from  small-scale  volumetric  inhomogeneities,  a  ray’s  interaction 
with  an  interface  across  which  seismic  velocities  and/or  material  densities  abruptly 
change  can  be  thought  of  as  a  form  of  scattering  radiation  pattern.  For  a  given  incident 
ray,  there  are  up  to  four  different  outgoing  ray  directions,  and  two  different  ray  types  (P 
and  S),  with  special  care  being  necessary  to  correctly  transform  the  S-polarization  angle. 
Aki  and  Richards  (1980)  develop  a  set  of  20  scattering  coefficients  that  cover  all 
possibilities  between  incoming  P  and  S  waves  and  their  associated  transmitted  or 
reflected  P  or  S  outgoing  waves.  These  coefficients  quantify  the  amplitude  and  phase  of 
the  outgoing  waves,  which  is  used  by  Radiative3D  to  compute  probabilities  of  an 
incoming  ray  transforming  into  the  corresponding  outgoing  ray  type.  This  allows  for  the 
modeling  of  reflection  (at  a  free  surface  or  interface)  and  transmission  (interface)  of  rays 
with  corresponding  possibility  of  transformation  between  P  and  S  polarization  types. 
This  feature  has  been  recently  completed  in  the  Radiative3D  codebase. 
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3.4.6  Features  in  Development 

Software  development  is  ongoing,  and  while  we  have  reached  a  state  where 
certain  studies  of  real-world  significance  can  be  undertaken  with  the  code  as-is,  there  still 
remain  some  features  being  tested  before  we  reach  full  target  functionality.  The  most 
significant  of  these  are: 

Tetrahedral  Model  Cells: 

Three-dimensional  Earth  models,  having  fully  general  spatial  gradients  of  velocity 
and  density  within  layers  need  to  be  discretized  in  some  fashion.  The  method  we  have 
chosen  is  to  divide  the  model  into  tetrahedral  cells  inside  of  which  the  Earth  parameters 
are  analytically  specified.  A  tetrahedron  is  the  polyhedron  with  the  minimum  number  of 
vertices  needed  to  specifiy  a  three-dimensional  volume,  which  allows  for  straighforward 
approaches  to  dividing  space  into  cellular  components.  Another  advantage  of  tetrahedra, 
(over,  e.g.,  rectangular  prisms),  is  that  Earth  parameters  specified  on  the  four  vertices 
allows  unique  and  complete  specification  of  the  properties  in  terms  of  linear  gradients 
throughout  the  volume  of  the  cell. 

Much  of  the  infrastructure  for  cellular  model  discretization  is  already  complete  in 
the  code,  and  is  utilized  in  the  layer-celled  structure  already  supported  in  the  code. 
Implementing  cells  with  tetrahedral  geometry  remains  primarily  an  implementation  detail 
at  this  point. 

Velocity  Gradients: 

In  a  uniform-velocity  medium,  ray  paths  are  straight  lines.  In  the  real  Earth,  rays 
can  follow  curvilinear  paths.  The  plan  for  Radiative3D  is  for  velocities  to  be  specified  as 
first-order  gradients,  resulting  in  paths  which  are  circular  arcs  withing  a  given  model  cell. 
This  allows  for  a  much  more  natural  representation  of  real  Earth  structure,  and  lessens  the 
detail  needed  to  represent  that  structure.  At  present  Radiative3D  supports  only  uniform 
velocities  (gradients  on  the  macro  level  can  be  represented  as  a  multicellular  stair-step 
model).  The  next  development  step  of  Radiative3D  will  be  to  implement  support  for 
velocity  gradients.  This  step  will  naturally  follow  the  implementation  of  tetrahedral 
model  cells,  since  this  geometry  allows  the  unique  specification  of  the  gradients  by 
specifying  known  velocities  on  the  cell  vertices. 

Anisotropy  of  Heterogeneity  Scale  Lengths 

Anisotropy  of  heterogeneity  scale  lengths  that  scatter  high  frequency  seismic 
waves  can  affect  coda  shapes  (Hong  and  Wu,  2005)  and  spatial  coherence  of  amplitudes 
(Tkalcic  et  al.,  2010)  depending  on  the  angle  their  wavefronts  make  with  the  axes  of 
vertically  or  horizontally  stretched  heterogeneity.  Nielsen  and  Thybo  (2006),  for 
example,  found  the  need  to  incorporate  horizontally  stretched  heterogeneity  to  model  Pn 
and  Sn  codas  observed  in  Russian  PNE  data.  We  have  made  some  tests  for  the  effects  of 


Approved  for  public  release;  distribution  is  unlimited. 

13 


the  radiation  patterns  of  vertically  or  horizontally  stretched  heterogeneity  (Cormier  et  al., 
201 1).  A  more  efficient  simulation  of  the  effects  of  heterogeneity  stretched  in  a  quasi¬ 
horizontal  direction  might  be  to  simply  combine  random  perturbations  to  thin  layering 
together  with  isotropic  heterogeneity.  This  can  be  accomplished  with  our  existing  code. 


3.4.7  Features  Not  Currently  Targeted: 

In  addition  to  the  code  features  detailed  above,  another  feature,  not  currently 
targetted  for  development,  remains  within  feasible  reach  of  the  codebase  as  it  now  stands, 
and  may  provide  a  starting  point  for  future  work.  In  particular,  while  Radiative3D 
currently  produces  seismic  envelope  plots,  it  is  technically  possible  to  produce  full- 
waveform  seismograms  by  tracking  not  only  the  amplitude,  but  also  the  phase,  of  each 
phonon  as  it  propagates  through  the  model.  By  binning  the  seismic  energy  as  amplitude 
phasors  rather  than  energy  scalars,  full  seismograms,  rather  than  envelopes,  can  be 
produced.  Phase  tracking  through  direct  propagation  is  trivial,  as  it  is  proportional  to  the 
product  of  frequency  and  travel  time.  Slightly  more  intricate  are  the  phase  shifts  that 
result  from  scattering  events,  and  from  reflection/transmission  interactions  with 
discontinuity  or  free-surface  interfaces.  These  phase  factors,  however,  are  actually 
already  calculated  by  the  codebase  when  it  calculates  scattering  and 
reflection/transmission  probabilities.  Adding  the  accounting  necessary  to  track 
cumulative  phase  through  a  phonon  lifetime  would  be  a  very  tractable  problem. 

4.  RESULTS  AND  DISCUSSION 

4.1  Lop  Nor  Region:  Deterministic  Model  and  Data 

We  focused  our  modeling  efforts  on  the  Lop  Nor,  China  region  (Figure  6)  because  of 
it  is  one  of  few  regions  in  which  nuclear  tests  and  earthquake  waveforms  are  well 
recorded  at  accessible  global  seismic  stations  and  networks.  Sykes  and  Nettles  (2009) 
found  that  more  than  half  of  the  earthquakes  in  the  Reviewed  Event  Bulletin  (REB)  of  the 
IMS  that  occurred  within  100  km  of  six  Lop  Nor  test  sites  from  2000  through  2008. 
LANL  has  also  constructed  maps  of  regional  wave  propagation  efficiencies  in  this  region 
(Phillips  et  al.,  1999).  Lop  Nor  is  located  near  the  southeastern  side  of  the  Tien  Shan,  a 
region  of  moderate  earthquake  activity  and  contemporary  horizontal  compressive  stress 
in  the  earth’s  crust,  which  can  be  an  important  factors  in  inducing  tectonic  release  from 
underground  nuclear  tests  in  this  region. 
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Figure  6.  Lop  Nor  test  site  region  and  paths  to  seismic  stations  MAK  and  WUS.  Lop  Nor 
test  site  and  paths  to  regional  seismic  stations  MAK  and  WUS  on  which  we  have 
concentrated  modeling  tests. 


The  Chinese  station  Urumqi  (assigned  the  station  code  WMQ  by  seismologists)  is 
about  2.15°  (250  km)  from  Lop  Nor.  Stations  MAK  and  WUS  are  ~  6.85°  from  Lop  Nor. 
For  the  earthquake  and  explosion  data  we  use  seismograms  from  events  downloaded 
from  IRIS-DMC. 

Figure  7  shows  some  waveform  examples  from  the  Lop  Nor  region.  The  data 
come  from  an  explosion  at  the  Lop  Nor  site  as  recorded  at  station  MAK  at  great  circle 
distance  6.85°  and  an  earthquake  in  the  same  area  recorded  at  both  MAK  and  WUS. 

WUS  is  approximately  at  the  same  distance  from  Lop  Nor  as  MAK  but  in  at  different 
azimuth.  The  left  column  in  Figure  7  shows  the  data  bandpassed  between  1-2  Hz  and  the 
right  column  shows  a  bandpass  between  6-8  Hz.  Not  only  are  there  differences  between 
the  earthquake  and  explosion  data  at  MAK,  there  are  also  differences  between  the 
earthquake  traces  at  MAK  and  WUS  in  the  time  windows  commonly  used  to  measure 
Pn/Lg  and  Pg/Lg  amplitude  ratios. 
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1  -2  Hz.  Explosion, Lop  Nor-MAK  6-8  Hz,  Explosion, Lop  Nor-MAK 


1  -2  Hz,  Earthquake,  Lop  Nor-MAK 


6-8  Hz,  Earthquake,  Lop  Nor-MAK 


1  -2  Hz,  Earthquake,  Lop  Nor-WUS  6-8  Hz,  Earthquake,  Lop  Nor-WUS 


Figure  7.  Vertical  component  bandpassed  data  recorded  at  WUS  and  MAK  from 
the  Lop  Nor  region. 


Using  known  elevations  of  Lop  Nor  and  seismic  stations  MAK  and  WUS,  along  with 
Moho  depths  from  the  Cornell  Moho  model  (atlas.geo.comell.edu/geoid/imagegrid.htmli 
and  layer  profiles  from  CRUST  2.0  (Laske  et  al.,  201 1)  at  those  same  locations,  we 
located  and  oriented  a  set  of  five  crust  layers  (sediments,  upper,  middle,  and  lower  crust, 
and  top  layer  of  the  mantle).  Additionally,  we  defined  an  additional  16  mantle  layers 
from  AK-135-F  (Kennett  et  ah,  1995)  and  regionalized  perturbations  (Gudmundson  and 
Sambridge,  1998)  to  a  depth  of  859  km.  These  layers  served  as  a  model  in  which  to  ran 
early  test  synthetics. 
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Figure  8.  Simplified  Lop  Nor  region  earth  model.  A  thin  sediment  layer  is  included. 


4.2  Sensitivity  of  Coda  Envelopes  to  Statistical  Parameters 

Since  there  are  no  documented,  large-scale,  differences  in  crustal  thicknesses 
between  the  paths  to  MAK  and  WUS,  and  we  could  not  find  any  previous  work  that  has 
connected  observed  coda  with  parameters  describing  a  heterogeneity  spectrum,  we  have 
performed  a  series  of  experiments  to  examine  the  sensitivity  of  regional  phases,  in 
general,  to  parameters  that  describe  the  small-scale  statistical  structure.  All  tests  of  the 
effects  of  statistical  parameters  have  employed  the  von  Karman  power  spectral  density 
function  (PSDF)  (Figure  1).  The  advantage  of  this  spectrum  is  that  it  can  incorporate 
changes  in  slope  of  the  power  spectrum  over  a  broad  band  of  wavenumbers,  typical  of 
media  having  a  fractal-like  behavior  of  scale  lengths,  as  well  as  a  comer  wavenumber 
above  which  the  power  spectrum  more  strongly  decays,  consistent  with  either  the 
existence  of  a  limiting  small  scale  or  a  limiting  sensitivity  of  the  wavefield  at  frequencies 
outside  the  pass  band  of  a  filtered  seismogram.  The  parameters  varied  are  the  fractional 

fluctuation  s  of  wave  velocities,  (£  =  =  —  ),  and  a  parameter  v  controlling  the 

Vp  Vs 

density  perturbation  as  factor  multiplying  the  P  velocity  fluctuation,  a  scale  length  a  that 
controls  the  wavenumber  comer,  and  a  von  Karman  order  number  a:  that  controls  the 
slope  of  the  power  spectrum  after  its  comer. 

We  investigated  the  sensitivity  of  regional  phase  codas  to  heterogeneity 
parameters  using  two  methods:  (a)  calculation  of  the  mean  free  path  (MFP)  and  the 
dipole  projection  (DP)  while  varying  each  of  the  parameters  and  holding  the  others 
constant  and  (b)  simulations  of  earthquakes  and  explosions  at  a  number  of  frequencies 
and  a  variety  of  parameters. 
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The  MFP  and  DP  calculations  provided  a  fairly  quick  and  accurate  way  to 
understand  the  general  effects  of  the  parameters.  The  MFP  reports  how  far  a  phonon  will 
travel  (in  km)  before  it  encounters  a  scattering  event.  Longer  MFPs  mean  less  scattering. 
The  dipole  projection  gives  information  about  the  direction  the  phonon  takes  after  it  is 
scattered  and  varies  from  -1  to  1 .  A  value  of  1  means  it  is  100%  likely  to  continue  in  its 
forward  direction,  and  a  value  of  -1  means  it  is  100%  likely  to  be  totally  scattered 
backwards. 

In  figures  plotting  the  effects  of  varying  a  parameter,  MFP  is  plotted  on  the  left 
side  y-axis  and  shown  as  solid  lines,  while  DP  is  plotted  on  the  right  side  y-axis  and 
shown  as  dashed  lines,  and  the  parameter  being  varied  is  plotted  on  the  x-axis.  The  red 
lines  show  the  effect  of  scattering  on  the  P-wave  phonons,  and  the  blue  lines  show  S- 
wave  phonons.  For  all  parameter  variations,  synthetic  envelopes  are  computed  at  2  Hz. 
The  effects  of  varying  each  heterogeneity  parameter  on  the  synthetic  envelopes  for 
earthquakes  and  explosions  are  compared  to  the  effects  of  a  reference  set  of  parameters 
(Table  1).  Qp  is  assumed  to  be  greater  than  Qs  by  a  factor  of  approximately  9/4, 
consistent  with  an  assumption  of  viscoelastic  attenuation  purely  in  shear  and  a  Poisson's 
ratio  of  %  (e.g.,  Anderson,  1989).  Figure  9  shows  the  synthetic  envelope  for  the 
reference  parameters  listed  in  Table  1  with  some  of  the  features  labeled.  For  comparison 
with  our  mantle  values,  Hedlin,  et  al.  (1997)  report  the  P  velocity  perturbation 
throughout  the  mantle  as  approximately  1%  and  for  a  equal  to  8  km.  The  frequency 
band  of  their  teleseismic  data,  however,  is  not  sensitive  to  the  smaller  scale  length  a 
required  needed  to  model  the  coda  of  regional  higher  frequency  data. 


Table  1.  Reference  heterogeneity  parameters  used  for  synthetic  comparisons. 
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Figure  9.  Coda  envelope  synthetics  for  earthquakes  in  the  Lop  Nor  region  for  a 
reference  heterogeneity  spectrum.  The  vertical  dotted  red  line  gives  the  location  in  time 
of  the  energy  maximum,  which  is  also  listed  at  the  top  of  the  figure.  The  bottom  left 
corner  of  the  figure  includes  the  frequency  of  the  run  (in  this  case  2  Hz)  and  the  isotropy 
of  the  moment  tensor  representation  of  the  source  (0%  indicates  earthquakes  and  100% 
indicates  explosions).  Phonon  statistics  listed  on  the  right  side  of  the  figure  include  the 
total  number  of  phonons  cast  during  this  simulation  (in  this  case  140  million),  the  total 
phonons  caught  at  this  location,  and  the  catch  rate.  Also  listed  on  the  right  side  is  the 
sum  of  the  energy  captured.  The  total  energy  is  indicated  as  P+S. 

4.2.1  Magnitude  of  Velocity  Perturbation  8 

Figure  10  shows  the  effects  of  velocity  perturbations  on  MFP  and  DP.  The 
dashed  lines  reflect  a  constant  S-wave  dipole  projection  of  0.48  and  a  constant  P-wave 
dipole  projection  of  0.18.  Varying  epsilon  has  no  effect  on  the  scattered  phonon’s 
preferred  direction  of  travel  after  scattering  because  the  magnitude  of  epsilon  does  not 
affect  the  shape  of  the  scattering  radiation  pattern  at  a  scattering  event.  On  the  other  hand, 
there  is  an  apparent  exponential  decrease  in  the  MFP  for  both  P  and  S  phonons  for  this 
range  of  velocity  perturbations.  We  should  expect  to  see  increased  coda  production  for 
both  S  and  P  waves  in  the  synthetic  envelopes  with  increased  epsilon. 
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Figure  10.  Mean  free  path  and  dipole  projection  for  velocity  perturbations 

Figures  1  la,b  compare  a  high  velocity  variation  (epsilon)  in  the  crust  to  the 
Reference  synthetic  envelopes  (0.08  vs.  0.04)  for  a  frequency  of  2  Hz.  For  the 
earthquake,  the  maximum  energy  and  maximum  phonon  count  in  the  Reference  ar  e  both 
associated  with  the  Lg  phase  (occurring  at  ~  241  s).  For  the  explosion,  the  maximum 
energy  and  maximum  phonon  count  in  the  Reference  are  both  associated  with  the  P  coda 
(occiming  at  ~  130  s).  We  see  a  number  of  effects  in  the  seismic  envelopes  caused  by  the 
change  in  the  heterogeneity  parameter,  including  change  in  total  energy7  and  phonons 
arriving  at  the  receiver,  change  in  maximum  energy  (giving  an  indication  of  phase 
amplitude),  onset  of  P  and  Lg  phases,  shape  of  direct  phase  arrival  (pulse),  and  coda 
following  the  direct  arrival. 
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Figure  1 1  a.  Earthquake  (left)  and  explosion  (right)  coda  envelopes  for  low  epsilon. 
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Figure  1  lb.  Earthquake  (left)  and  explosion  (right)  coda  envelopes  for  high  epsilon. 


The  combined  results  of  the  MFP/DP  plots  and  synthetic  envelopes  show  that 
with  higher  velocity  perturbations,  the  phonons  have  smaller  mean  free  paths,  thus 
encountering  more  scattering  events.  The  lack  of  change  in  the  dipole  polarization, 
however,  demonstrates  that  energy  tends  to  scatter  more  towards  the  direction  of  the 
paths  before  the  scattering  events,  and  will  eventually  reach  their  original  destination. 

The  shift  in  arrival  time  of  the  maximum  energy  in  the  earthquake  envelope  from  the  Lg 
arrival  window  to  the  P  arrival  window  reflects  the  fact  that  S-waves  are  more  sensitive 
to  velocity  fluctuations  than  P  waves. 

The  scattering  caused  by  the  increase  in  the  velocity  perturbation  produces  a 
marked  redistribution  in  energy,  stretching  it  out  over  a  longer  period  of  time.  This  is 
represented  by  a  slower  decay  of  coda  envelopes,  the  broadening  of  the  pulses,  and  a 
pronounced  delay  in  the  phonon  arrivals  (shaded  area)  These  results  reflect  the  increased 
production  of  scattered  P  and  S  energy  as  well  as  conversions  between  and  P  and  S 
energy  (particularly  for  the  explosion).  Without  scattering,  almost  all  of  the  explosion 
energy  is  in  the  P  arrival.  With  the  scattering,  in  addition  to  coda  energy  following  the  P- 
arrival,  there  is  significant  energy  arriving  at  the  time  we  would  expect  to  see  Lg.  This 
delayed  energy,  however,  is  reflected  in  a  slower  coda  envelope  decay  following  the  P 
phases,  which  cannot  easily  be  mistaken  for  the  typical  Lg  coda  envelope  observed  from 
earthquakes.  This  strongly  suggests  that  mechanisms  other  than  scattering  alone  need  to 
be  considered  to  explain  observations  of  anomalous  low  Pg/Lg  ratios  from  some 
explosions. 

4.2.2  Magnitude  of  Density  Perturbation  v 


The  density  perturbation  is  a  multiplier  of  the  velocity  perturbation:—  =  v^-. 

Figure  12  shows  the  effects  of  density  perturbations  from  0.3  to  1.5  on  MFP  and  DP. 
Higher  values  of  nu  produce  a  shorter  mean  free  path  as  expected.  In  contrast  to  the 
effect  of  increasing  velocity  variations,  the  effect  of  increasing  density  variations  is  a 
linear  rather  than  the  exponential  decrease  in  mean  free  path.  The  higher  values  of  nu 
also  lower  the  dipole  projection,  which  means  the  phonons  are  scattered  further  away 
from  the  forward  direction  of  travel.  This  suggests  that  nu  variations  will  be  important  in 


Approved  for  public  release;  distribution  is  unlimited. 
21 


controlling  tlie  attenuation  of  the  peak  amplitudes  of  regional  phases  defined  from  narrow 
windows  of  group  velocity. 
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Figure  12.  Mean  free  path  and  dipole  projection  for  density  perturbations. 


Figure  13  compares  a  high  density  perturbation  to  the  Reference  synthetic 
envelopes  (1.6  vs.  0.8  in  the  crust).  The  maximum  energy  occurs  at  the  same  time  as  the 
Reference,  and  the  maximum  phonon  count  occurs  only  slightly  later  (~  255  s).  The  P 
onset  appears  to  be  sharper.  There  is  a  lower  amplitude  as  some  of  the  energy  arrival  is 
delayed,  but  the  effect  is  not  as  pronounced  as  that  of  the  velocity  perturbation. 
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Figure  13  Synthetic  coda  envelopes  for  density  perturbations.  Top:  Reference  earthquake 
(left)  and  explosion  (right).  Bottom:  effects  of  increased  density  perturbation. 
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In  general,  the  influence  of  density  variations  is  to  produce  more  scattered  energy 
than  for  velocity  variations  alone.  The  ability  to  adjust  nu  could  be  important  in  other 
seismic  coda  studies,  such  as  those  involving  partial  melt  (Hong,  et  al.,  2004).  Shearer 
and  Earle  (2008)  take  nu  0.8  for  small-scale  heterogeneity  (1-10  km).  For  larger  scale 
structure  in  the  mantle,  geodynamicists  take  nu  on  the  order  of  0.2  or  less  because  large 
values  cannot  be  sustained  by  buoyancy. 


4.2.3  Scale  length  a 


Figure  14  shows  the  effects  of  scale  length  variations  from  0.5  to  12  km  on  MFP 
and  DP.  For  the  velocities  and  frequency  (2  Hz)  used  in  coda  envelope  simulations,  the 
wavelength  in  the  crust  is  2  km  for  the  S  waves  and  3.4  km  for  P  waves. 
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Figure  14.  Mean  fee  path  and  dipole  projection  for  variable  scale  lengths  of  a  von 
Kaiman  spectrum. 

Figure  15  shows  synthetic  envelopes  for  three  changes  in  the  scale  length: 
doubled  compared  to  the  Reference  (0.4  km),  3  km,  and  12  km.  This  allowed  us  to 
examine  the  effects  of  the  scale  length  being  greater  than  the  wavelength.  The  mean  free 
path  exponentially  decreases  as  the  scale  length  approaches  the  wavelength,  but  the 
dipole  projection  exponentially  increases  toward  1  (forward  scattering).  Not  surprisingly, 
the  competition  between  these  two  effects  results  in  the  changes  the  shape  and 
attenuation  of  regional  phase  codas  most  strongly  when  the  scale  length  is  on  the  order  of 
the  wavelength,  near  the  intersection  of  the  MFP  and  DP  curves  in  Figure  13. 
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Location:  799.39  km  at  318.68  deg  az  from  source 
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Figure  15.  Synthetic  coda  envelopes  for  variable  scale  lengths  of  a  von  Karman 
spectrum,  (a.)  earthquake  Reference,  b.)  explosion  Reference,  c.)  earthquake  0.4  km  scale 
length,  d.)  explosion  0.4  km  scale  length,  e.)  earthquake  3  km  scale  length,  f.)  explosion  3 
km  scale  length,  g.)  earthquake  12  km  scale  length,  h.)  explosion  12  km  scale  length). 
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4.2.4  Hurst  Parameter  k 

The  fluctuation  spectra  for  the  von  Kaiman  medium  are  fiat  up  to  a  comer  wave 
number  inversely  proportional  to  the  con  elation  distance  and  then  fall  off  at  higher  wave 
numbers  (ka  >  1).  It  is  also  characterized  by  heterogeneities  that  are  self-similar  for  ka  > 

1 .  This  medium  is  "rougher"  at  small  length  scales  than  the  exponential  medium  (Frankel 
and  Clayton,  1986). 

Figure  16  shows  the  effects  of  k  variations  from  0. 1  to  0.8  on  MFP  and  DP.  As 
kappa  increases,  the  MFP  decreases,  which  cause  scattering  delays  in  the  energy  arrival. 
As  the  slope  of  the  heterogeneity  spectrum  flattens  above  the  comer  wavenumber  the 
scattering  is  increasingly  in  the  forward  direction  even  though  scattering  events  are  more 
frequent.  The  effect  of  increased  scattering  at  higher  kappa  wifi  have  less  effect  on 
redistributing  energy  and  stretching  coda  envelopes  but  may  have  a  stronger  effect  on 
equilibrating  energy  on  the  three  components  of  motion  due  to  a  more  fr  equent  sampling 
of  the  heterogeneity  radiation  patterns  that  control  polarization  and  wave  type  (P  or  S). 
The  heterogeneity  power,  however,  will  decrease  faster  at  wavenumbers  above  the  comer 
wavenumber  for  higher  kappa.  Thus  the  effect  of  different  k’s  can  be  expected  to  be 
minor,  which  we  found  to  be  true  in  our  experiments  with  synthetic  coda.  Since  visible 
effects  are  not  easily  discernible,  we  do  not  show  synthetics  for  this  experiment. 
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Figure  16.  Mean  free  path  and  dipole  projection  for  the  Hurst  parameter  k,  which  controls 
spectral  shape  for  wavelengths  larger  than  the  scale  length  of  a  von  Kaiman  spectrum. 
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4.3  Intrinsic  Q  Effects 


Radiative3D  allows  us  to  treat  intrinsic,  viscoelastic  attenuation  separately  from 
scattering  to  help  determine  their  relative  contributions  to  the  total  attenuation  of  peak 
coda  amplitudes  or  to  separate  their  contributions  to  measured  coda-Q’s.  Figure  17 
compares  coda  envelopes  against  the  synthetic  envelopes  for  the  Reference  deterministic 
and  statistical  model  of  heterogeneity  with  and  without  intrinsic  attenuation.  For  the 
earthquake,  the  2  Hz.  total  energy  without  Q"1  included  is  about  100  times  that  of  the 
reference  envelope  with  the  maximum  energy  value  about  50  times  that  of  the  reference. 
The  phonon  capture  is  about  the  same.  The  effect  of  Q'1  is  to  reduce  the  amount  of  energy 
in  each  phonon  as  a  function  of  time  traveled.  For  the  explosion,  the  2  Hz  total 
energy, without  Q"1  included,  is  about  10  times  that  of  the  reference  synthetics,  with  the 
maximum  energy  value  about  twice  that  of  the  reference.  The  phonon  capture  is  about  the 
same. 
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Figure  17.  Synthetic  coda  envelopes  for  earthquakes  and  explosions  with  and  without 
intrinsic  attenuation.  Top  left:  earthquake  reference;  top  right:  explosion  reference; 
bottom  left:  earthquake  without  intrinsic  attenuation;  bottom  right  explosion  without 
intrinsic  attenuation. 
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This  test  shows  that  not  adding  intrinsic  attenuation  results  in  the  coda  energy  being 
unrealistically  extended  in  time.  Future  work  with  the  radiative  transport  algorithm  will 
allow  the  parameter  space  of  scattering  intrinsic  attenuation  to  be  fully  explored  to 
examine  tradeoffs  that  can  produce  a  fit  to  regional  coda  shapes.  Accurate  models  of 
intrinsic  Q  in  the  crust  and  just  beneath  the  Moho  are  essential  for  proper  modeling  of  Sn 
and  Lg. 

4.4  Surface  layer  effects 

A  highly  heterogeneous,  strongly  scattering,  sediment  layer  may  affect  regional 
phase  coda  (Baumgardt,  2001),  particularly  for  shallow  emplaced  explosions.  Figure  18 
shows  example  synthetic  envelopes  exploring  the  effects  of  a  0.5  km  thick  sediments 
layer  with  the  velocities,  density,  and  Q’s  given  in  Table  2.  The  first  row  of  the  envelopes 
are  the  Reference  synthetics  from  the  heterogeneity  parameter  study,  hi  the  second  row  of 
envelopes  the  velocities  in  the  sediments  layer  have  been  adjusted  to  nearly  match  those 
of  the  crust  layer  immediately  below.  The  thin  layer  in  which  low  velocities  have  been 
inactivated  is  still  treated  separately  for  scattering,  although  for  this  test  is  epsilon  equal 
to  0.04,  matching  the  values  used  in  the  thicker  crustal  layer. 


Table  2.  Deterministic  properties  of  test  sediment  layer. 


VP 

Vs 

Density 

Qp 

Qs 

Sedi  Layer 

2.50 

1.20 

2.10 

162.8 

50.0 

Upper  Crust  Layer 

6.13 

3.53 

2.75 

2261.7 

1000.0 

For  both  the  earthquake  and  the  explosion,  the  most  obvious  result  of  the  low 
velocity  sediment  layer  is  on  the  P/Lg  amplitude  ratio.  Although  the  differences  between 
earthquake  and  explosion  envelopes  are  strong  in  both  cases,  the  effect  of  the  sediment 
layer  will  be  to  slightly  reduce  the  effectiveness  of  the  P/Lg  discriminant.  The 
instantaneous  state  of  scattered  phonons  arriving  at  the  receiver  will  be  more  equal  in  P 
and  S  for  the  case  of  the  sediment  layer.  This  is  also  apparent  in  the  tendency  of  the 
sediment  scattering  layer  to  equilibrate  energy  on  all  three  components  of  motion  (X,  Y, 
Z). 
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Figure  18.  Synthetic  coda  envelopes  for  earthquakes  and  explosions  with  and  without  a 
thin  sediment  layer.  Top  left:  earthquake  Reference;  top  right:  explosion  Reference; 
bottom  left:  earthquake  no  sediments;  bottom  right:  earthquake  no  sediments. 

4.5  Pn  and  Sn 

High  frequency  Pn  and  Sn  are  traditionally  thought  of  as  head  waves  traveling 
just  beneath  the  Moho  discontinuity.  Early  attempts  to  model  high  frequency  Pn  and  Sn 
as  classical  head  waves  failed  because  the  amplitude  in  a  classical  head  decays  with 
frequency  as  1/omega.  These  phases  are  now  recognized  as  interference  head  waves 
(Cerveny  and  Ravindia.  1971),  which  have  a  representation  as  series  of  body  waves  that 
are  multiply  reflected  along  the  underside  for  either  a  Moho  following  Earth’s  curvature 
or  for  velocities  that  increase  with  depth  beneath  the  Moho.  Our  initial  experiments 
assumed  a  flat  earth  and  homogeneous  layers,  separated  by  at  most  planar  tilted  layers.  In 
these  models  energy  could  not  return  from  the  mantle  except  as  scattered  waves.  With  a 
series  of  stair  step  discontinuities,  we  have  simulated  the  effects  of  Earth  curvature  via  an 
earth-flattening  tr  ansformation  and  positive  velocity  gradients  with  decreasing  radius 
(Figure  19).  These  modifications  allow  Pn  and  Sn  to  return  to  the  surface  after  being 
multiply  reflected  by  and  scattered  beneath  the  Moho.  Envelopes  were  synthesized  for 
the  southern  Xinjiang  earthquake  recorded  at  MAK  (Figure  20)  for  this  Moho  model. 
Figure  2 1  shows  the  results  of  this  experiment  with  scattering  from  the  reference 
statistical  model  turned  on  and  off.  A  sharp  Pn  arrival  is  visible  in  Figure  2 1 .  Sharp 
arrivals  preceding  Lg,  which  can  be  interpreted  as  Sn,  are  best  visible  in  the  envelopes 
having  scattering  turned  off. 
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Cross  Section  from  LOP  to  MAK 
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Figure  19.  Cross-section  of  an  Earth  model  having  a  Moho  transition  specified  by  a 
series  of  stair-step  layers. 


Figure  20.  Lop  Nor  paths  to  MAK  and  WUS  and  source  mechanism  of  a  southern 
Xinjiang  earthquake.  Left:  receiver  placement,  phonon  gather  radii,  and  epicenter  of 
southern  Xinjiang  earthquake.  Right:  mechanism  of  earthquake. 
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Location:  799.39  km  at  318.68  deg  az  from  source 
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Figure  21.  Synthetic  coda  envelopes  for  earthquakes  and  explosions  in  a  Lop  Nor  model 
having  a  Moho  transition.  Top:  Synthetic  envelopes  for  the  Moho  transition  model  for 
the  path  to  station  MAK  from  the  southern  Xianjiang  earthquake  mechanism  shown  in 
Figure  22  with  scattering  predicted  from  the  reference  statistical  model  turned  on. 
Bottom:  synthetic  envelopes  for  the  same  path  with  scattering  turned  off 


Figure  22  compares  synthetic  envelope  against  the  recorded  earthquake.  There  is 
good  agreement  in  arrival  time  and  relative  amplitude  for  the  Pn,  Pg,  and  Lg  phases.  The 
Sn  phase  is  much  smaller  in  the  synthetic  compared  to  the  data.  We  believe  that  reducing 
the  intrinsic  S  attenuation  may  allow  more  of  the  Sn  energy  to  show  up  for  this  phase. 


Station  MAK,  Z  channel  @  2.0  Hz,  2003-03-13  15:07:06 


Figure  22.  Comparison  of  2  Hz  synthetics  with  observed  waveforms  and  passed  for  2  Hz.  Comparison 
of  data  (below)  with  synthetic  coda  envelope  (above)  for  a  model  having  a  Moho  transition  layer. 
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4.6  Synthetics  compared  with  data  by  frequency  band 

Radiative3D  is  mil  at  a  single  frequency  and  can  create  simulated 
seismogram  envelopes  at  any  frequency.  This  allows  analysts  to  interpret  the 
structural  or  source  origin  of  frequency- specifics  effects  on  narrow  bandpassed 
seismic  tr  aces.  Table  3  lists  the  values  of  the  parameters  used  to  generate  synthetic 
envelopes  that  include  Pn  and  Sn  phases.  These  are  similar  to  the  Reference 
parameters  used  in  the  heterogeneity  parameter  sensitivity  study  with  a  transition 
layer  added  at  the  Moho  as  discussed  above. 


Table  3.  Scattering  par  ameters  for  fr  equency  com 
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eps 
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kappa 
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Crust 
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300 

There  is  now  a  substantial  body  of  empirical  evidence  that  regional  P/S  ratios 
provide  poor  discrimination  below  some  fr  equency,  typically  about  2  to  3  Hz,  and  useful 
discrimination  at  higher  frequencies,  e.g.,  Fisk  2006.  Hence,  in  Figure  23  we  compare 
computed  envelopes  against  narrow  bandpassed  data  at  2,  3,  and  4  Hz. 


Figure  23.  Synthetic  coda  envelopes  compared  with  data  recorded  at  MAK  for  a  series 
of  frequencies.  Synthetics  compared  against  data  recorded  at  station  MAK  for  the  2003- 
03-13  earthquake  for  frequencies  2,  3,  and  4  Hz.  The  approximate  arrivals  of  the 
regional  crustal  phases  (Pg  and  Lg)  are  marked.  The  relative  changes  in  amplitude  of 
the  Pg  and  Lg  phases  produced  by  the  synthetics  agree  quite  well  with  the  data. 
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4.7  Tectonic  release:  regional  propagation  effects 

We  initiated  an  investigation  of  the  effect  on  high  frequency  regional  phase 
propagation  of  explosion  sources  associated  with  components  of  tectonic  stress  release. 
These  experiments  were  inspired  in  part  by  the  work  of  Bukchin  et  al.,  2001.  That  study 
jointly  inverted  P  wave  first  motions  and  teleseismic  Love  and  Rayleigh  wave  spectra  for 
unconstrained  moment  tensors  of  earthquakes  and  nuclear  tests  in  the  Lop  Nor  region. 
The  amount  of  tectonic  release  of  nuclear  tests  was  expressed  as  an  angle  given  by  the 
inverse  tangent  of  the  ratio  of  the  isotropic  to  non-isotropic  moment  tensor  components. 
The  experimental  series  incorporated  moment  tensors  reported  by  Bukchin  et  al.  for  the 
Lop  Nor  1996-06-08  nuclear  tests,  although  the  source  location  used  was  that  of  the 
southern  Xinjian  2003-03-13  earthquake.  Envelopes  were  synthesized  for  paths  to 
stations  MAK  and  WUS.  Additional  event  mechanisms  and  locations,  including 
experiments  in  which  depths  of  tectonic  release  vary,  are  available  on  our  group’s  wiki 
web  link  (rainbow.phys.uconn.edu/geowiki). 

For  this  series  of  experiments,  we  used  the  moment  tensor  decompositions  of 
several  explosion  events  published  by  Bukchin  et  al.  as  our  starting  points.  Although  they 
treat  the  isotropic  and  deviatoric  components  as  vertically  separated,  in  these  runs  we 
treated  them  as  collocated  and  coherent  (meaning  a  single  moment  tensor  specifies  the 
entire  source,  as  opposed  to  running  the  deviatoric  and  isotropic  source  radiation  events 
separately  and  summing  the  synthetics.  We  found  that  separations  of  the  explosion  and 
tectonic  release  components  as  large  of  4  km  did  not  visibly  affect  coda  envelopes. 

Figure  24  shows  an  example  of  one  of  our  event  simulations  compared  to  data 
recorded  in  several  narrow  bandpassed  filters  for  1  to  4  Hz.  Note  that  the  P/Lg  ratio  for 
the  explosion  with  tectonic  release  is  largest  in  the  lowest  frequency  band  of  1  Hz.  This 
ratio,  however,  is  also  typical  of  that  found  for  explosions  without  tectonic  release  at  1 
Hz,  even  though  this  example  is  one  of  the  nuclear  tests  for  which  Bukchin  et  al. 
determined  the  largest  amount  of  tectonic  release,  primarily  constrained  from  teleseismic 
surface  waves.  In  higher  frequency  bands,  excitation  of  P  coda  is  much  larger  than  Lg, 
suggesting  that  event  discrimination  based  on  P/Lg  ratios  will  not  suffer  in  bands  2  Hz 
and  higher.  This  behavior  will  be  strongly  affected  by  the  intrinsic  Q  for  shear  waves  in 
the  crust.  We  have  assumed  a  QP=  2000  and  Qs  =  1000  in  the  crust.  Even  though  the 
assumed  shear  Q  is  quite  high  it  is  still  low  enough  that  any  Lg  trapped  in  the  crust  from 
the  explosion  is  strongly  attenuated  at  this  distance  range  at  frequencies  2  Hz  and  higher. 
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Figure  24.  Synthetic  coda  envelopes  for  nuclear  tests  with  tectonic  release  compared  with 
observed  data  for  Lop  Nor  nuclear  tests  at  a  series  of  narrow  frequency  bands.  Synthetic 
envelopes  (left)  and  observed  data  (right)  for  Lop  Nor  nuclear  tests  as  a  function  of  frequency. 
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5.  CONCLUSIONS 


We  have  developed  a  radiative  transport  code  to  shoot  body  wave  rays  through 
general  deterministic  3-D  structure,  including  the  computation  of  quantities  required  to 
synthesize  high  frequency  body  wave  coda  generated  by  small-scale,  statistically 
described  heterogeneity.  Representing  the  wavefield  as  a  sum  of  multiply  scattered  P  and 
S  waves  in  3-D,  it  includes  reflection/transmission/P-S  conversion  by  interfaces,  effects 
on  the  amplitude  and  polarization  of  scattered  re -radiated  waves  by  from  statistically 
described  small-scale  heterogeneity,  and  intrinsic  viscoelastic  attenuation.  Velocity 
gradients  and  earth  curvature  are  handled  by  the  ability  to  include  thin  layers  and  an 
earth- flattening  transformation.  Development  of  an  option  to  handle  complex  3-D 
gradients  and  generalized  boundary  topography  in  the  large-scale  background  model  is 
nearly  complete,  with  velocity  gradients  and  topography  defined  from  model  values  at 
tetrahedral  knots. 


Coda  envelopes  can  be  plotted  in  several  styles:  (1)  as  a  function  of  time  and 
component  of  motion  at  a  single  station,  (2)  as  a  function  of  distance  and  time  to  assist 
in  travel-time  picking  and  travel-time  uncertainty  estimates,  or  (3)  as  snapshots  in  time 
as  a  function  of  depth  and  range  to  understand  the  evolution  of  P  and  S  energy  and  the 
homogenization  of  radiation  pattern. 


Source  radiation  patterns  can  be  input  in  the  form  of  generalized  moment  tensors, 
which  can  include  tectonic  release.  The  effect  of  source  spectra  shape  on  coda  envelopes 
can  be  handled  by  varying  the  scalar  moment  for  each  frequency  in  simulations  for 
individual  narrow  bandpassed  simulation. 

The  availability  of  both  explosion  and  earthquake  waveforms,  detailed  maps  of 
regional  phase  efficiency,  and  the  existence  of  3-D  deterministic  structural  models 
derived  from  global  stations  and  transportable  arrays  make  the  Lop  Nor  region  an  ideal 
region  to  implement  our  modeling  technique.  In  this  region  we  have  completed 
modeling  experiments  that  have  tested  the  effects  on  regional  coda  of  parameters 
specifying  the  heterogeneity  spectrum,  the  effects  of  crustal  and  mantle  intrinsic 
attenuation,  the  effect  of  a  thin  sediment  layer,  the  nature  of  the  Moho  transition,  and  the 
effects  and  detection  of  explosion  triggered  tectonic  release.  Tests  varying  the 
parameters  specifying  small-scale  heterogeneity  find  that  the  strongest  effects  on  regional 
codas  occur  when  the  scale  length  specifying  the  comer  of  a  von  Karman  spectrum  are 
close  to  the  dominant  wavelength  of  a  narrow  bandpassed  simulation.  Pn  and  Sn 
amplitudes  strongly  depend  on  velocity  gradients  and  structural  complexity  beneath  the 
Moho.  Accurate  estimates  of  intrinsic  Q  in  the  crust  are  critical  input  to  the  modeling  of 
Lg  and  the  prediction  of  P/Lg  ratios.  Tectonic  release  detected  from  low  frequency, 
teleseismic,  surface  waves  minimally  affects  the  high  frequency  (>2  Hz)  regional  coda  in 
the  Lop  Nor  region.  Thus  far,  the  most  significant  conclusion  of  all  of  these  tests  is  that 
performance  of  discriminants  based  on  P/Lg  amplitude  ratios  is  best  at  frequencies  above 
2  Hz  from  recordings  at  ranges  less  than  1000  km. 
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